With your quantum computers and scanners now working the way they are supposed to, you are ready to identify the optimal planet for your gravity assist maneuver. A well-executed slingshot may just be the key to your escape.
Your optimization protocol reveals one clear winner: a purple gas giant, with glimmering gaseous swirls and a large spot that reminds you of Jupiter. This planet is much different from what your broken systems would’ve identified if you hadn’t repaired them!
As you prepare your ship for the slingshot, you discover a major obstacle: the region around the planet is rife with space debris. Myriads of objects of varying sizes, some as small as sand grains, some as large as beach balls. Metals and silicates and rocks, remnants of asteroid collisions and perhaps even civilizations long-gone. The objects revolve in an expansive region around the planet.
Your slingshot is at risk: with the debris orbiting at speeds over twenty-five thousand kilometers per hour, a collision could be disastrous.
But, you do not lose hope. You have quantum computers that work, you know quantum optimization, and you have a drone on board that could effectively collect any trash you target.
Complete the exercises below in order to clear an optimal path through the space debris by successfully collecting specific orbiting objects with your drone, thus mitigating the risk associated with the slingshot.
To clear our pathway around the planet, we will be sending our drone out to the debris to collect the major pieces; the pieces that would cause catastrophic damage if our starship collided with them. With our scanners now in tip-top shape, we are able to estimate the distances between these large pieces of debris. We cannot risk any time wasted; the black hole is pulling us closer and closer. We need to calculate the optimal pathway for our drone to take in order to collect all pieces of debris in the smallest amount of time (corresponding to the shortest total path for this mission).
The collection of the debris using the shortest path is essentially a well-known problem called the traveling salesman problem. In addition to being a notorious NP-complete problem that has drawn the attention of computer scientists and mathematicians for over two centuries, the Traveling Salesman Problem (TSP) has important bearings on finance and marketing, as its name suggests. Colloquially speaking, the traveling salesman is a person that goes from city to city to sell merchandise. The objective in this case is to find the shortest path that would enable the salesman to visit all the cities and return to their hometown, i.e. the city where they started traveling. By doing this, the salesman gets to maximize potential sales in the least amount of time.
This problem derives its importance from its "hardness" and ubiquitous equivalence to other relevant combinatorial optimization problems that arise in practice.
The mathematical formulation with some early analysis was proposed by W.R. Hamilton in the early 19th century. Mathematically the problem is, as in the case of Max-Cut [1], best abstracted in terms of graphs. The TSP on the nodes of a graph asks for the shortest Hamiltonian cycle that can be taken through each of the nodes. A Hamilton cycle is a closed path that uses every vertex of a graph once. The general solution is unknown and an algorithm that finds it efficiently (e.g., in polynomial time) is not expected to exist.
Find the shortest Hamiltonian cycle in a graph $G=(V,E)$ with $n=|V|$ nodes and distances, $w_{ij}$ (distance from vertex $i$ to vertex $j$). A Hamiltonian cycle is described by $N^2$ variables $x_{i,p}$, where $i$ represents the node and $p$ represents its order in a prospective cycle. The decision variable takes the value 1 if the solution occurs at node $i$ at time order $p$. We require that every node can only appear once in the cycle, and for each time a node has to occur. This amounts to the two constraints (here and in the following, whenever not specified, the summands run over 0,1,...N-1)
$$ \begin{equation}\tag{1} \sum_{i} x_{i,p} = 1 ~~\forall p. \end{equation} $$
$$ \begin{equation}\tag{2} \sum_{p} x_{i,p} = 1 ~~\forall i. \end{equation} $$
For nodes in our prospective ordering, if $x_{i,p}$ and $x_{j,p+1}$ are both 1, then there should be an energy penalty if $(i,j) \notin E$ (not connected in the graph). The form of this penalty is
$$\sum_{i,j\notin E}\sum_{p} x_{i,p}x_{j,p+1}>0,$$where it is assumed the boundary condition of the Hamiltonian cycles $(p=N)\equiv (p=0)$. However, here it will be assumed a fully connected graph and not include this term. The distance that needs to be minimized is
$$C(\textbf{x})=\sum_{i,j}w_{ij}\sum_{p} x_{i,p}x_{j,p+1}.$$Putting this all together in a single objective function to be minimized, we get the following:
$$C(\textbf{x})=\sum_{i,j}w_{ij}\sum_{p} x_{i,p}x_{j,p+1}+ A\sum_p\left(1- \sum_i x_{i,p}\right)^2+A\sum_i\left(1- \sum_p x_{i,p}\right)^2,$$where $A$ is a free parameter. One needs to ensure that $A$ is large enough so that these constraints are respected. One way to do this is to choose $A$ such that $A > \mathrm{max}(w_{ij})$.
Once again, it is easy to map the problem in this form to a quantum computer, and the solution will be found by minimizing an Ising Hamiltonian. [2]
from itertools import permutations
import matplotlib.axes as axes
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
from typing import Optional
# import all necessary qiskit packages
from qiskit.algorithms import MinimumEigensolver, NumPyMinimumEigensolver, VQE
from qiskit.algorithms.minimum_eigensolvers import QAOA, SamplingVQE
from qiskit.algorithms.optimizers import L_BFGS_B, SLSQP, SPSA
from qiskit.circuit.library import TwoLocal
from qiskit.tools.visualization import plot_histogram
from qiskit.utils import algorithm_globals, QuantumInstance
from qiskit_aer import Aer
from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit_optimization.applications import Maxcut, Tsp
from qiskit_optimization.problems import QuadraticProgram
<frozen importlib._bootstrap>:219: RuntimeWarning: scipy._lib.messagestream.MessageStream size changed, may indicate binary incompatibility. Expected 56 from C header, got 64 from PyObject
# Draw and visualize graph solution
def draw_graph(G, colors, pos):
default_axes = plt.axes(frameon=True)
nx.draw_networkx(G, node_color=colors, node_size=900, alpha=0.9, ax=default_axes, pos=pos, node_shape="o")
edge_labels = nx.get_edge_attributes(G, "weight")
nx.draw_networkx_edge_labels(G, pos=pos, edge_labels=edge_labels)
# Compute cost of the obtained result and display result
def TSPCost(energy, result_bitstring, adj_matrix):
print("-------------------")
print("Energy:", energy)
print("Tsp objective:", energy + offset)
print("Feasibility:", qubo.is_feasible(result_bitstring))
solution_vector = tsp.interpret(result_bitstring)
print("Solution Vector:", solution_vector)
print("Solution Objective:", tsp.tsp_value(solution_vector, adj_matrix))
print("-------------------")
# Draw Graph
draw_tsp_solution(tsp.graph, solution_vector, colors, pos)
# Plots a convergence graph
def PlotGraph(ideal, mean):
plt.figure(figsize=(8, 6))
# Plot label
for value in range(len(mean)):
plt.plot(mean[value], label="Proposal {}".format(value+1) )
# Ideal plot
plt.axhline(y=ideal, color="tab:red", ls="--", label="Target")
plt.legend(loc="best")
plt.xlabel("Optimizer iteration")
plt.ylabel("Energy")
# Plot graph title
plt.title("TSP Entangler Line {} Constraint".format(len(mean)))
plt.show()
# Let us first generate a graph with 4 nodes!
#### enter your code below ####
n = 4
num_qubits = 16
#### enter your code above ####
tsp = Tsp.create_random_instance(n, seed=123)
adj_matrix = nx.to_numpy_matrix(tsp.graph)
print("distance\n", adj_matrix)
colors = ["r" for node in tsp.graph.nodes]
pos = [tsp.graph.nodes[node]["pos"] for node in tsp.graph.nodes]
draw_graph(tsp.graph, colors, pos)
distance [[ 0. 48. 91. 33.] [48. 0. 63. 71.] [91. 63. 0. 92.] [33. 71. 92. 0.]]
from qc_grader.challenges.fall_2022 import grade_lab3_ex1
grade_lab3_ex1(n, num_qubits) # Expected result type: int, int
Submitting your answer. Please wait... Congratulations 🎉! Your answer is correct and has been submitted.
Great start! You have successfully created a graph that shows the debris location and their distances from one another.
The brute force approach involves trying every combination through classical computation to find the optimal solution. As you can infer, this procedure becomes more computationally-intensive as the number of nodes increases, or as the number of locations the traveling salesman must visit, increases.
from itertools import permutations
def brute_force_tsp(w, N):
a = list(permutations(range(1, N)))
last_best_distance = 1e10
for i in a:
distance = 0
pre_j = 0
for j in i:
distance = distance + w[j, pre_j]
pre_j = j
distance = distance + w[pre_j, 0]
order = (0,) + i
if distance < last_best_distance:
best_order = order
last_best_distance = distance
print("order = " + str(order) + " Distance = " + str(distance))
return last_best_distance, best_order
best_distance, best_order = brute_force_tsp(adj_matrix, n)
print(
"Best order from brute force = "
+ str(best_order)
+ " with total distance = "
+ str(best_distance)
)
def draw_tsp_solution(G, order, colors, pos):
G2 = nx.DiGraph()
G2.add_nodes_from(G)
n = len(order)
for i in range(n):
j = (i + 1) % n
G2.add_edge(order[i], order[j], weight=G[order[i]][order[j]]["weight"])
default_axes = plt.axes(frameon=True)
nx.draw_networkx(
G2, node_color=colors, edge_color="b", node_size=600, alpha=0.8, ax=default_axes, pos=pos
)
edge_labels = nx.get_edge_attributes(G2, "weight")
nx.draw_networkx_edge_labels(G2, pos, font_color="b", edge_labels=edge_labels)
draw_tsp_solution(tsp.graph, best_order, colors, pos)
order = (0, 1, 2, 3) Distance = 236.0 Best order from brute force = (0, 1, 2, 3) with total distance = 236.0
To solve this optimization problem, we first need to convert the TSP problem into a quadratic program, a program which lays out the optimization problem as a quadratic objective function and relevant constraint functions. The TSP defines variables as binary variables, so you will see those defined in the quadratic program as well.
qp = tsp.to_quadratic_program()
print(qp.prettyprint())
Problem name: TSP
Minimize
48*x_0_0*x_1_1 + 48*x_0_0*x_1_3 + 91*x_0_0*x_2_1 + 91*x_0_0*x_2_3
+ 33*x_0_0*x_3_1 + 33*x_0_0*x_3_3 + 48*x_0_1*x_1_0 + 48*x_0_1*x_1_2
+ 91*x_0_1*x_2_0 + 91*x_0_1*x_2_2 + 33*x_0_1*x_3_0 + 33*x_0_1*x_3_2
+ 48*x_0_2*x_1_1 + 48*x_0_2*x_1_3 + 91*x_0_2*x_2_1 + 91*x_0_2*x_2_3
+ 33*x_0_2*x_3_1 + 33*x_0_2*x_3_3 + 48*x_0_3*x_1_0 + 48*x_0_3*x_1_2
+ 91*x_0_3*x_2_0 + 91*x_0_3*x_2_2 + 33*x_0_3*x_3_0 + 33*x_0_3*x_3_2
+ 63*x_1_0*x_2_1 + 63*x_1_0*x_2_3 + 71*x_1_0*x_3_1 + 71*x_1_0*x_3_3
+ 63*x_1_1*x_2_0 + 63*x_1_1*x_2_2 + 71*x_1_1*x_3_0 + 71*x_1_1*x_3_2
+ 63*x_1_2*x_2_1 + 63*x_1_2*x_2_3 + 71*x_1_2*x_3_1 + 71*x_1_2*x_3_3
+ 63*x_1_3*x_2_0 + 63*x_1_3*x_2_2 + 71*x_1_3*x_3_0 + 71*x_1_3*x_3_2
+ 92*x_2_0*x_3_1 + 92*x_2_0*x_3_3 + 92*x_2_1*x_3_0 + 92*x_2_1*x_3_2
+ 92*x_2_2*x_3_1 + 92*x_2_2*x_3_3 + 92*x_2_3*x_3_0 + 92*x_2_3*x_3_2
Subject to
Linear constraints (8)
x_0_0 + x_0_1 + x_0_2 + x_0_3 == 1 'c0'
x_1_0 + x_1_1 + x_1_2 + x_1_3 == 1 'c1'
x_2_0 + x_2_1 + x_2_2 + x_2_3 == 1 'c2'
x_3_0 + x_3_1 + x_3_2 + x_3_3 == 1 'c3'
x_0_0 + x_1_0 + x_2_0 + x_3_0 == 1 'c4'
x_0_1 + x_1_1 + x_2_1 + x_3_1 == 1 'c5'
x_0_2 + x_1_2 + x_2_2 + x_3_2 == 1 'c6'
x_0_3 + x_1_3 + x_2_3 + x_3_3 == 1 'c7'
Binary variables (16)
x_0_0 x_0_1 x_0_2 x_0_3 x_1_0 x_1_1 x_1_2 x_1_3 x_2_0 x_2_1 x_2_2 x_2_3
x_3_0 x_3_1 x_3_2 x_3_3
To solve this problem, we need to convert our quadratic program into a form that a quantum computer can work with. In this case we will generate an Ising Hamiltonian, which is simply a representation (i.e. model) of the energy of this particular system (i.e. problem). We can use eigensolvers to find the minimum energy of the Ising hamiltonian, which corresponds to the shortest path and the solution of our problem.
How do we get the Ising Hamiltonian? First, we convert the quadratic program to a QUBO (Quadratic Unconstrained Binary Optimization), a combinatorial optimization problem representation to then obtain an Ising Hamiltonian. The next exercise will teach you how we do this.
from qiskit_optimization.converters import QuadraticProgramToQubo
#### enter your code below ####
qp2qubo = QuadraticProgramToQubo()
qubo = qp2qubo.convert(qp)
qubitOp, offset = qubo.to_ising()
#### enter your code above ####
print("Offset:", offset)
print("Ising Hamiltonian:", str(qubitOp))
Offset: 51756.0 Ising Hamiltonian: -6456.0 * IIIIIIIIIIIIIIIZ - 6456.0 * IIIIIIIIIIIIIIZI - 6456.0 * IIIIIIIIIIIIIZII - 6456.0 * IIIIIIIIIIIIZIII - 6461.0 * IIIIIIIIIIIZIIII - 6461.0 * IIIIIIIIIIZIIIII - 6461.0 * IIIIIIIIIZIIIIII - 6461.0 * IIIIIIIIZIIIIIII - 6493.0 * IIIIIIIZIIIIIIII - 6493.0 * IIIIIIZIIIIIIIII - 6493.0 * IIIIIZIIIIIIIIII - 6493.0 * IIIIZIIIIIIIIIII - 6468.0 * IIIZIIIIIIIIIIII - 6468.0 * IIZIIIIIIIIIIIII - 6468.0 * IZIIIIIIIIIIIIII - 6468.0 * ZIIIIIIIIIIIIIII + 1592.5 * IIIIIIIIIIIIIIZZ + 1592.5 * IIIIIIIIIIIIIZIZ + 1592.5 * IIIIIIIIIIIIIZZI + 1592.5 * IIIIIIIIIIIIZIIZ + 1592.5 * IIIIIIIIIIIIZIZI + 1592.5 * IIIIIIIIIIIIZZII + 1592.5 * IIIIIIIIIIIZIIIZ + 12.0 * IIIIIIIIIIIZIIZI + 12.0 * IIIIIIIIIIIZZIII + 12.0 * IIIIIIIIIIZIIIIZ + 1592.5 * IIIIIIIIIIZIIIZI + 12.0 * IIIIIIIIIIZIIZII + 1592.5 * IIIIIIIIIIZZIIII + 12.0 * IIIIIIIIIZIIIIZI + 1592.5 * IIIIIIIIIZIIIZII + 12.0 * IIIIIIIIIZIIZIII + 1592.5 * IIIIIIIIIZIZIIII + 1592.5 * IIIIIIIIIZZIIIII + 12.0 * IIIIIIIIZIIIIIIZ + 12.0 * IIIIIIIIZIIIIZII + 1592.5 * IIIIIIIIZIIIZIII + 1592.5 * IIIIIIIIZIIZIIII + 1592.5 * IIIIIIIIZIZIIIII + 1592.5 * IIIIIIIIZZIIIIII + 1592.5 * IIIIIIIZIIIIIIIZ + 22.75 * IIIIIIIZIIIIIIZI + 22.75 * IIIIIIIZIIIIZIII + 1592.5 * IIIIIIIZIIIZIIII + 15.75 * IIIIIIIZIIZIIIII + 15.75 * IIIIIIIZZIIIIIII + 22.75 * IIIIIIZIIIIIIIIZ + 1592.5 * IIIIIIZIIIIIIIZI + 22.75 * IIIIIIZIIIIIIZII + 15.75 * IIIIIIZIIIIZIIII + 1592.5 * IIIIIIZIIIZIIIII + 15.75 * IIIIIIZIIZIIIIII + 1592.5 * IIIIIIZZIIIIIIII + 22.75 * IIIIIZIIIIIIIIZI + 1592.5 * IIIIIZIIIIIIIZII + 22.75 * IIIIIZIIIIIIZIII + 15.75 * IIIIIZIIIIZIIIII + 1592.5 * IIIIIZIIIZIIIIII + 15.75 * IIIIIZIIZIIIIIII + 1592.5 * IIIIIZIZIIIIIIII + 1592.5 * IIIIIZZIIIIIIIII + 22.75 * IIIIZIIIIIIIIIIZ + 22.75 * IIIIZIIIIIIIIZII + 1592.5 * IIIIZIIIIIIIZIII + 15.75 * IIIIZIIIIIIZIIII + 15.75 * IIIIZIIIIZIIIIII + 1592.5 * IIIIZIIIZIIIIIII + 1592.5 * IIIIZIIZIIIIIIII + 1592.5 * IIIIZIZIIIIIIIII + 1592.5 * IIIIZZIIIIIIIIII + 1592.5 * IIIZIIIIIIIIIIIZ + 8.25 * IIIZIIIIIIIIIIZI + 8.25 * IIIZIIIIIIIIZIII + 1592.5 * IIIZIIIIIIIZIIII + 17.75 * IIIZIIIIIIZIIIII + 17.75 * IIIZIIIIZIIIIIII + 1592.5 * IIIZIIIZIIIIIIII + 23.0 * IIIZIIZIIIIIIIII + 23.0 * IIIZZIIIIIIIIIII + 8.25 * IIZIIIIIIIIIIIIZ + 1592.5 * IIZIIIIIIIIIIIZI + 8.25 * IIZIIIIIIIIIIZII + 17.75 * IIZIIIIIIIIZIIII + 1592.5 * IIZIIIIIIIZIIIII + 17.75 * IIZIIIIIIZIIIIII + 23.0 * IIZIIIIZIIIIIIII + 1592.5 * IIZIIIZIIIIIIIII + 23.0 * IIZIIZIIIIIIIIII + 1592.5 * IIZZIIIIIIIIIIII + 8.25 * IZIIIIIIIIIIIIZI + 1592.5 * IZIIIIIIIIIIIZII + 8.25 * IZIIIIIIIIIIZIII + 17.75 * IZIIIIIIIIZIIIII + 1592.5 * IZIIIIIIIZIIIIII + 17.75 * IZIIIIIIZIIIIIII + 23.0 * IZIIIIZIIIIIIIII + 1592.5 * IZIIIZIIIIIIIIII + 23.0 * IZIIZIIIIIIIIIII + 1592.5 * IZIZIIIIIIIIIIII + 1592.5 * IZZIIIIIIIIIIIII + 8.25 * ZIIIIIIIIIIIIIIZ + 8.25 * ZIIIIIIIIIIIIZII + 1592.5 * ZIIIIIIIIIIIZIII + 17.75 * ZIIIIIIIIIIZIIII + 17.75 * ZIIIIIIIIZIIIIII + 1592.5 * ZIIIIIIIZIIIIIII + 23.0 * ZIIIIIIZIIIIIIII + 23.0 * ZIIIIZIIIIIIIIII + 1592.5 * ZIIIZIIIIIIIIIII + 1592.5 * ZIIZIIIIIIIIIIII + 1592.5 * ZIZIIIIIIIIIIIII + 1592.5 * ZZIIIIIIIIIIIIII
from qc_grader.challenges.fall_2022 import grade_lab3_ex2
grade_lab3_ex2(qubitOp) # Expected result type: PauliSumOp object
Submitting your answer. Please wait... Congratulations 🎉! Your answer is correct and has been submitted.
Good one!
Let's first attempt to obtain our reference solution by solving our QUBO with an exact classical method, NumPyMinimumEigensolver:
# Making the Hamiltonian in its full form and getting the lowest eigenvalue and eigenvector
ee = NumPyMinimumEigensolver()
result = ee.compute_minimum_eigenvalue(qubitOp)
energy_numpy = result.eigenvalue.real
# Call helper function to compute cost of the obtained result and display result
TSPCost(energy = result.eigenvalue.real, result_bitstring = tsp.sample_most_likely(result.eigenstate), adj_matrix = adj_matrix)
------------------- Energy: -51520.0 Tsp objective: 236.0 Feasibility: True Solution Vector: [1, 2, 3, 0] Solution Objective: 236.0 -------------------
Here we are now going to try running the same problem on a quantum computer.
However, since the calculations for this particular problem can take some time (10 minutes or even longer), we recommend that you take a look at one of the solutions your crew member obtained earlier when running the same problem. We know that time is precious. Especially in a scenario where we have a blackhole to escape from.
This is not a graded exercise so you can take a look at the results and proceed to the next segment of this lab.
Here is one of the results your crew member obtained when running the same problem earlier. The results can be slightly different from when you run the same problem due to the heuristic nature of the algorithm, but you can see that it is not showing good convergence to the solution. The 'solution objective' (solution of our objective function) gave us '302.0,' which is higher than our reference solution (shortest path) '236.0.'
Note: Since the cell below will take some time finish running, we recommend that you run this later on your spare time. Luckily, one of your crew members has already run the same problem on a quantum computer earlier and saved the results for you so that you don't have to execute. You may proceed to the next segment.
# # #This is not a graded exercise. Run this cell when you have spsare time.
# # # Let us first specify the backend we will be running
# ### Specifying the backend and creating quantum instance
# algorithm_globals.random_seed = 123
# seed = 10598
# backend = Aer.get_backend("qasm_simulator")
# quantum_instance = QuantumInstance(backend, seed_simulator=seed, seed_transpiler=seed)
# ### Building an ansatz with rotational Y gates
# # We can use a traditional VQE with an ansatz, which is a trial function selected as an initial guess approximating the minimum eigenstate of a Hamiltonia, built with Y single-qubit rotations and entangler steps.
# # Please note that this is heuristic ansatz, which means that it provides you with an educated guess to help solve a problem but does not guarantee to be the most optimal answer.
# spsa = SPSA(maxiter=300)
# ry = TwoLocal(qubitOp.num_qubits, "ry", "cz", reps=5, entanglement="linear")
# vqe = VQE(ry, optimizer=spsa, quantum_instance=quantum_instance)
# result = vqe.compute_minimum_eigenvalue(qubitOp)
# ### Print solution
# z = tsp.interpret(x)
# print("solution:", z)
# print("solution objective:", tsp.tsp_value(z, adj_matrix))
# draw_tsp_solution(tsp.graph, z, colors, pos)
In many applications it is important to find the minimum eigenvalue of a matrix. For example, in chemistry, the minimum eigenvalue of a Hermitian matrix characterizing the molecule is the ground state energy of that system. As we learned earlier, we can map our shortest path problem into an Ising Hamiltonian, which allows us to search for the shortest path in the same way we find the ground state energy of a molecule.
You may have heard that a quantum state can be represented as a vector in a complex vector space. This space where all the possible quantum states exist is mathematically called a 'Hilbert Space.' As your optimization problem grows in size and the more number of constraints you need to take into consideration, this adds a degree of complexity that will grow your Hilbert Space exponentially.
Problems like the Traveling Salesman Problem are known to be NP-Hard due to how our solution space grows exponentially with each city and constraints added.
This 'NP-hardness' in theoretical computer science make heuristics the only viable option for a variety of complex optimization problems that need to be routinely solved in real-world applications where the objective of a heuristic approach is to produce a solution in a reasonable time frame that is good enough for solving the problem at hand. This solution may not be the best of all the solutions to this problem, or it may simply approximate the exact solution. But it is still valuable because finding it does not require a prohibitively long time.
You may realize that this indicates that there is a 'trade-off' between optimality and execution time. Some heuristics converge faster than others. How can we get better convergence with our current quantum systems with limited number of qubits? In Part II, we will explore tnis notion by looking at an approach that can help reduce our search space so we can arrive at a good solution in a relatively reasonable amount of time.
The “natural form” of a QUBO model inherently is an unconstrained problem, other than those requiring the variables to be binary (hence the name Quadratic Unconstrained Binary Optimization). Most of the problems that are of interest to mankind tend to include additional constraints in the formulation of the problem to obtain good solutions. Many of these constrained models can be effectively re-formulated as a QUBO model by introducing penalties into the objective function as an alternative to explicitly imposing constraints to best fit the problem construction. [3]
We can utilize the constraints of our TSP problem to dynamically construct a problem-specific PQC that reflects those constraints of the optimization problem. Therefore, we can restrict a unitary transformation that is provided by the problem-specific PQC while taking constraints into account. This makes it possible to make the search space smaller in hopes of arriving at an optimal or near-optimal solution!
As mentioned above, we shall be looking at implementing a problem specific PQC for our TSP problem to attempt to reduce our search space and arrive at a good solution. This is an active area of interest and there are multiple PQC approaches for the TSP problem, however specifically, we shall be looking at one of the implementations by Matsuo, A., Suzuki, Y., & Yamashita, S. Problem-specific parameterized quantum circuits of the VQE algorithm for optimization problems. 2020. [4]
One of the reasons of intractability for many problems comes from the huge search space that exists for a solution. With so many possibilities to evaluate, a brute force search approach quickly becomes computationally unfeasible. The idea here is to reduce the subspace in such a way that it still includes the solution subspace, but is significantly smaller than the complete subspace, lets call it $|S_{all}|$. So the main objective is to propose a PQC that gives us a subspace: $S_{proposed}$ that includes $S_{feasible\_solutions}$, but such that the size of $S_{proposed}$ subspace is smaller than $|S_{all}|$. This increases the possibility of arriving at optimal or near optimal solution while reducing the required time to find the same as compared to, let's say, using a heuristic approach on the whole subspace!
Usually, an optimization problem can have more than one constraint. For such cases, we create multiple problem-specific parameterized quantum sub-circuits, each of which reflect the corresponding constraint. Then, by combining those sub-circuits properly, even though the optimization problem has more than one constraint, it is still possible to create a problem-specific PQC and reduce the search space.
The paper we mentioned above proposes 4 different types of PQCs approaching the problem in different ways, each with different characteristics and constraints considered. For the NISQ era, characteristics of a circuit such as number of gates, parameters, entanglement and circuit depth are important factors when designing PQCs and here we shall be looking into a few of the approaches as showcased in the paper above. We will be utilizing the VQE algorithm to build upon the routine to get our results using the Estimator and Sampler primitive!
Let's check out the first type of problem-specific PQC as proposed for our TSP problem here to effectively collect our space debris. The paper proposes us to take into account only the constraints in the first line as shown below:
$$\sum_{i}^{N} x_{i,p} = 1 ~~\forall p = 1...N.$$As we can see, in each of the constraints, exactly one of the variables has to be one while the other variables have to be zero. This type of constraint restricts the set of solutions to the set of bases corresponding to the $W$ state.
A $W$ state is a superposition of states where exactly one of the qubits is $|1\rangle$ while all other qubits are $|0\rangle$ with equal amplitudes. A $W$ state for $n$ qubits is represented as $|W\rangle = \frac{1}{\sqrt{2^n}} (|10...0\rangle + |01...0\rangle + |00...1\rangle)$. Each of the bases of this state corresponds to the constraint as mentioned above $\sum{x_i} = 1$. Since the other bases do not satisfy the constraint of $\sum{x_i} = 1$, we do not consider them here hence reducing our subspace for solving the problem with considerations of the constraints from the equation in the first line. We shall create circuits representing the $W$ state in such a way that we can control the amplitudes of each base.
$$|W(φ)\rangle = \sum_{i} α_{i(φ)}|ψ_i\rangle, $$$$\sum_{i} {|α_{i(φ)}|}^2 = 1, $$$$|ψ_i\rangle ∈ \{{|10...0\rangle, |01...0\rangle, |00...1\rangle}\}$$Let us see how we can do that for a 3 qubit example. Construct the W state for a 3 qubit circuit below!
# Define 3 Node graph for our learning exercises
n = 3
num_qubits = n**2
tsp = Tsp.create_random_instance(n, seed=250)
adj_matrix = nx.to_numpy_matrix(tsp.graph)
# Plot graph
colors = ["tab:red" for node in tsp.graph.nodes]
pos = [tsp.graph.nodes[node]["pos"] for node in tsp.graph.nodes]
print("Learning problem 3 Node graph: ")
draw_graph(tsp.graph, colors, pos)
# Define quadratic program for 4 node graph
qp = tsp.to_quadratic_program()
# Define Ising operator for 3 Node graph
qp2qubo = QuadraticProgramToQubo()
qubo = qp2qubo.convert(qp)
qubitOp, offset = qubo.to_ising()
# Making the Hamiltonian in its full form and getting the lowest eigenvalue and eigenvector
ee = NumPyMinimumEigensolver()
result = ee.compute_minimum_eigenvalue(qubitOp)
# Ideal result
ideal = result.eigenvalue.real
print(f"Ideal energy: {ideal}")
Learning problem 3 Node graph: Ideal energy: -4751.0
from qiskit import QuantumRegister, ClassicalRegister
from qiskit.circuit import ParameterVector, QuantumCircuit
# Initialize circuit
qr = QuantumRegister(3)
cr = ClassicalRegister(3)
circuit = QuantumCircuit(qr,cr)
# Initialize Parameter with θ
theta = ParameterVector('θ', 2)
#### enter your code below ####
circuit.x(0)
circuit.barrier()
circuit.ry(theta[0],1)
circuit.cz(0,1)
circuit.ry(-theta[0],1)
circuit.barrier()
circuit.ry(theta[1],2)
circuit.cz(1,2)
circuit.ry(-theta[1],2)
circuit.barrier()
circuit.cnot(1,0)
circuit.cnot(2,1)
# Apply gates to circuit
#### enter your code above ####
circuit.measure(qr,cr)
circuit.draw("mpl")
from qc_grader.challenges.fall_2022 import grade_lab3_ex3
grade_lab3_ex3(circuit) # Expected result type: QuantumCircuit object
Submitting your answer. Please wait... Congratulations 🎉! Your answer is correct and has been submitted.
Now its time to run the circuit with some parameters and observe the quasi-probability distribution! As per the paper we should be seeing the probability distribution being controlled as per the equation (11) in section 3.2.1 in the paper.
$$α_{1(φ)} |100\rangle + α_{2(φ)} |010\rangle + α_{3(φ)} |001\rangle,$$$$\sum_{i}^3 {|α_{i(φ)}|}^2 = 1, $$$$α_{1(φ)} = cos(\theta_1), α_{2(φ)} = -sin(\theta_1) cos(\theta_2), α_{3(φ)} = sin(\theta_1) sin(\theta_2)$$For an example's sake, let's try to get an equal distribution for the 3 qubit state. That would mean we should be seeing equal/near equal distributions for the states $|100\rangle, |010\rangle, |001\rangle$ i.e. decimal states $1, 2 $ and $4$ if we satisfy the equation above to have $α_{1(φ)} = α_{2(φ)} = α_{3(φ)}$
Let us visualize and confirm this result! According to the equation above, solving it to see a near equal probability across the three states of 1, 2 and 4, we can enter in $\theta_1$ as 54.73 degrees and $\theta_2$ as 45 degrees!
from qiskit_ibm_runtime import QiskitRuntimeService
# Initialize service and backend
#QiskitRuntimeService.save_account(channel='ibm_quantum', token='my_token', overwrite=True) #uncomment if you need to save your account again
service = QiskitRuntimeService(channel="ibm_quantum")
# Set simulator
backend = service.backends(simulator=True)[0]
print(backend)
<IBMBackend('ibmq_qasm_simulator')>
from math import pi
from qiskit_ibm_runtime import Session, Estimator, Sampler, Options
options = Options(simulator={"seed_simulator": 42},resilience_level=0) #DO NOT CHANGE
#### Enter your code below ####
# Define sampler object
with Session(service=service, backend=backend):
sampler = Sampler(options=options)
result = sampler.run(circuits=[circuit],parameter_values=[54.73*np.pi/180,45*np.pi/180],shots=5000).result()# Run the sampler. Remember, the theta here is in degrees! :)
# Plot result
result_dict = result.quasi_dists[0]# Obtain the quasi distribution. Note: Expected result type: QuasiDistribution dict
#### Enter your code above ####
values = list(result_dict.keys()) # Obtain all the values in a list
probabilities = list(result_dict.values()) # Obtain all the probabilities in a list
plt.bar(values,probabilities,tick_label=values)
plt.show()
print(result_dict)
{1: 0.3298, 4: 0.3286, 2: 0.3416}
from qc_grader.challenges.fall_2022 import grade_lab3_ex4
grade_lab3_ex4(result_dict) # Expected result type: QuasiDistribution dict
Submitting your answer. Please wait... Congratulations 🎉! Your answer is correct and has been submitted.
You may not see an exact distribution due to the randomness added by the simulator and the less precise theta value used here. But nonetheless, we have achieved our intended example objective! Feel free to play around with the theta values and as far as the thetas conform to the equation constraints we mentioned above, you should be seeing amplitude changes corresponding to only the states we mentioned above!
Now let's proceed to apply this to an algorithm and see how we can use this state for our problem as per the paper!
We will now proceed to build the routine as suggested above. For a TSP problem, the problem is encoded for an $N$ node graph is encoded in $N^2$ qubits. We need to append the W gate routine in this order for the PQC to be applied successfully as per our problem constraint. We will first look at a case for a 2 node graph and then proceed to have a generalized function to generate the W state for any arbitrary value of N.
For an example here, let us consider this condition for a 2 node graph. For a 2 node graph, we would require 4 qubits. Note, here we shall follow the row-major order. The first two qubits would represent the columns of the first row respectively as shown in the figure and the last 2 qubits would represent the columns in the $N^{th}$ row which in this case are the two elements in the second row. The number of parameters we require in general would be $(N)*(N-1)$. Let's check out how to apply this TSP entangler routine for a 2 node graph and extrapolate this for the general solution as an exercise.
# Two Node TSP entangler
from qiskit import QuantumRegister, ClassicalRegister
from qiskit.circuit import ParameterVector, QuantumCircuit
# Set value of N
N = 2
# Define Circuit and ParameterVector
circuit = QuantumCircuit(N**2)
theta = ParameterVector('theta',length=(N-1)*N)
# Block 1 ----------------
# X Gate -----------------
circuit.x(0)
circuit.barrier()
# Parameter Gates --------
circuit.ry(theta[0],1)
circuit.cz(0,1)
circuit.ry(-theta[0],1)
# CX gates ---------------
circuit.cx(1,0)
circuit.barrier()
# Block 2 ----------------
# X Gate -----------------
circuit.x(2)
circuit.barrier()
# Parameter Gates --------
circuit.ry(theta[1],3)
circuit.cz(2,3)
circuit.ry(-theta[1],3)
# CX gates ---------------
circuit.cx(3,2)
circuit.draw("mpl")
For a graph with $N = 3$, we will be following the same. The circuit you made in Exercise 3 will be the one that will be applied 3 times as a block. In general, you will have the W state applied for N qubits as one block. In total, there will be N blocks across $N^2$ qubits to complete the problem.
In the next exercise, we shall build up a TSP entangler routine for the above shown constraint for any general N qubit case. Pass the function to the grader to complete this exercise.
from qiskit import QuantumRegister, ClassicalRegister
from qiskit.circuit import ParameterVector, QuantumCircuit
def tsp_entangler(N):
# Define QuantumCircuit and ParameterVector size
circuit = QuantumCircuit(N**2)
theta = ParameterVector('theta',length=(N-1)*N)
#############
# Hint: One of the approaches that you can go for is having an outer `for` loop for each node and 2 inner loops for each block. /
# Continue on by applying the same routine as we did above for each node with inner `for` loops for /
# the ry-cz-ry routine and cx gates. Selecting which indices to apply for the routine will be the /
# key thing to solve here. The indices will depend on which parameter index you are checking for and the current node! /
# Reference solution has one outer loop for nodes and one inner loop each for rz-cz-ry and cx gates respectively.
#############
k=0
for i in range(0,N**2,N):
circuit.x(i)
for j in range(i+1,i+N):
circuit.ry(theta[k],j)
circuit.cz(j-1,j)
circuit.ry(-theta[k],j)
k+=1
for j in range(i,i+N-1):
circuit.cnot(j+1,j)
#### enter your code above ####
# X Gate
# Parameter Gates
# CX gates
#### enter your code above ####
return circuit
from qc_grader.challenges.fall_2022 import grade_lab3_ex5
grade_lab3_ex5(tsp_entangler) # Expected result type: function
Submitting your answer. Please wait... Congratulations 🎉! Your answer is correct and has been submitted.
Great job! We now have the building blocks to run our PQC and test our solution! Let's start constructing the VQE function for calculating our best result for the problem. We shall be using the VQE class to build up our VQE routine and calculate the result using the Estimator Primitive.
The VQE class in Qiskit is now refactored to leverage the Quantum primitive construct. It uses a variational technique to find the minimum eigenvalue of a given diagonal operator but is executed using the Estimator primitive.
An instance of VQE also requires an ansatz, a parameterized QuantumCircuit, to prepare the trial state $|\psi(\vec\theta)\rangle$. It also needs a classical optimizer which varies the circuit parameters $\vec\theta$ such that the expectation value of the operator on the corresponding state approaches a minimum. In this case we are interested in extracting the result bitstring corresponding to the optimal solution from the converged result.
These are the following values it takes as per the Qiskit Documentation:
Arguments:
estimator: The estimator primitive to compute the expectation value of the Hamiltonian operator.ansatz: A parameterized quantum circuit to prepare the trial state.optimizer: A classical optimizer to find the minimum energy. This can either be a Qiskit.Optimizer or a callable implementing the .Minimizer protocol.initial_point: An optional initial point (i.e. initial parameter values) for the optimizer. The length of the initial point must match the number of ansatz parameters. If None, a random point will be generated within certain parameter bounds. SamplingVQE will look to the ansatz for these bounds. If the ansatz does not specify bounds, bounds of -2\pi, 2\pi will be used.callback: A callback that can access the intermediate data at each optimization step. These data are: the evaluation count, the optimizer parameters for the ansatz, the estimated value, and the metadata dictionary.
Note on execution times:
The VQE cells may take up to 10-15 minutes to run based on your location. Here we will be passing in jobs to the simulator on cloud time here is due to each job taking about ~1 second to run including network latency of passing the jobs and the result back to the optimizer. Feel free to check out Program jobs on the IQX portal to see if your jobs are running: https://quantum-computing.ibm.com/jobs?jobs=runtime.
# Imports
from qiskit_optimization import QuadraticProgram
from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit.algorithms.minimum_eigen_solvers import NumPyMinimumEigensolver, MinimumEigensolver
from qiskit.algorithms.minimum_eigensolvers import SamplingVQE, VQE
from qiskit.circuit.library import TwoLocal
from qiskit.algorithms.optimizers import SPSA, SLSQP, L_BFGS_B
from qiskit.quantum_info import SparsePauliOp
from qiskit_ibm_runtime import Session, Estimator, Sampler, Options
import numpy as np
# Define PQC for 3 node graph (Hint: Call the previous function you just made)
PQC = tsp_entangler(3)#### enter your code here ####
def RunVQE(estimator, model, optimizer, operator, init=None):
# Store intermediate results
history = {"eval_count": [], "parameters": [], "mean": [], "metadata": []}
# Define callback function
def store_intermediate_result(eval_count, parameters, mean, metadata):
history["eval_count"].append(eval_count)
history["parameters"].append(parameters)
history["mean"].append(mean.real)
#### Enter your code below ####
# Define VQE run
vqe = VQE(estimator=estimator,ansatz=model,optimizer=optimizer,initial_point=init)
#
#
# build your code here
#
#
# Compute minimum_eigenvalue
result = vqe.compute_minimum_eigenvalue(operator)
#### Enter your code above ####
return result, history["mean"]
%%time
# Do not change the code below and use the given seeds for successful grading
from qiskit.utils import algorithm_globals
algorithm_globals.random_seed = 123
options = Options(simulator={"seed_simulator": 42},resilience_level=0)
# Define optimizer
optimizer = SPSA(maxiter=50)
# Define Initial point
np.random.seed(10)
init = np.random.rand((n-1)*n) * 2 * np.pi
# Define backend
backend = service.backends(simulator=True)[0]
#### Enter your code below ####
# Call RunVQE
with Session(service = service, backend = backend):
#----# Enter your code here. Do not forget to pass options when you pass the Estimator object
result, mean = RunVQE(estimator=Estimator(options=options),model=PQC,optimizer=optimizer,operator=qubitOp,init=init)
#### Enter your code above ####
# Print result
type(result)
CPU times: user 2.42 s, sys: 183 ms, total: 2.6 s Wall time: 6min 41s
qiskit.algorithms.minimum_eigensolvers.vqe.VQEResult
print(result)
{ 'aux_operators_evaluated': None,
'cost_function_evals': 100,
'eigenvalue': -4750.24525,
'optimal_circuit': <qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7fac1adf76a0>,
'optimal_parameters': { ParameterVectorElement(theta[1]): -1.5732974847551533,
ParameterVectorElement(theta[2]): 4.702581245089359,
ParameterVectorElement(theta[5]): 3.1488519394325616,
ParameterVectorElement(theta[4]): 4.712429783663164,
ParameterVectorElement(theta[0]): 3.1427472772537177,
ParameterVectorElement(theta[3]): 4.703140373187985},
'optimal_point': array([ 3.14274728, -1.57329748, 4.70258125, 4.70314037, 4.71242978,
3.14885194]),
'optimal_value': -4750.24525,
'optimizer_evals': None,
'optimizer_result': <qiskit.algorithms.optimizers.optimizer.OptimizerResult object at 0x7fac1ae5b280>,
'optimizer_time': 401.24259424209595}
from qc_grader.challenges.fall_2022 import grade_lab3_ex6
grade_lab3_ex6(result) # Expected result type: VQEResult object
Submitting your answer. Please wait... Congratulations 🎉! Your answer is correct and has been submitted.
Next here, we shall proceed to compute the Eigenstate of our solution and get our best distribution binary string from the VQEResult object. To obtain the eigenstate, we shall sample the parametrized ansatz with the optimal parameters that we obtained by running the Estimator call on our VQE function. We shall use the Sampler here to sample our ansatz with the optimal_parameters from the VQEResult object. The key corresponding to the maximum of the nearest_probability_distribution we obtain will be our binary state of the best solution which we will pass on to evaluate!
Let's setup our Sampler routine for the same. First, we will extract the optimal_circuit from our computed result and apply measurements for the Sampler routine. Here we use inplace=False to make sure we store the circuit with measurements as a new circuit.
optimal = result.optimal_circuit.measure_all(inplace=False) # Extract optimal circuit and add measurements
optimal.draw("mpl")
Next, let us proceed to setup our Sampler to obtain the best bitstring from the result we just computed. Here we shall pass the optimal circuit with the measurements as computed above and extract the optimal_parameters list from result.
A note on the usage of nearest_probability_distribution: While a Quasiprobability distribution does offer a wider gamut of information data points at your disposal, we cannot directly substitute it for a real probability distribution. For a Quasiprobability distribution to be converted into a true probability distribution, we shall use here the
nearest_probability_distribution()function. This function takes a quasiprobability distribution as an input and maps it to the closest probability distribution as defined by the L2-norm up to a certain error bound
We shall now use this to get our eigenstate and obtain the optimal result from this distribution!
Exercise 7:
Build up a Sampler routine and pass in the optimal circuit and the list of parameter values obtained from the VQE result we computed previously. Obtain the nearest probability distribution from the sampler result you just obtained. Please keep the seed values constant for the grading to be successful.
We shall use this value to compute the feasibility of our solution in the next cells.
def BestBitstring(result, optimal_circuit):
energy = result.eigenvalue.real
options = Options(simulator={"seed_simulator": 42},resilience_level=0) # Do not change and please pass it in the Sampler call for successful grading
#### enter your code below ####
with Session(service = service, backend = backend):
sampler_result = Sampler(options=options).run(circuits=optimal_circuit,parameter_values=list(result.optimal_parameters.values())).result()
result_prob_dist = sampler_result.quasi_dists[0].nearest_probability_distribution()
#### enter your code above ####
max_key = format(max(result_prob_dist, key = result_prob_dist.get),"016b")
result_bitstring = np.array(list(map(int, max_key)))
return energy, sampler_result, result_prob_dist, result_bitstring
%%time
# Compute values
with Session(service=service, backend=backend):
energy, sampler_result, result_prob_dist, result_bitstring = BestBitstring(result=result, optimal_circuit=optimal)
print("Optimal bitstring = ", result_bitstring[7:])
Optimal bitstring = [0 1 0 1 0 0 0 0 1] CPU times: user 51.9 ms, sys: 3.97 ms, total: 55.9 ms Wall time: 5.91 s
from qc_grader.challenges.fall_2022 import grade_lab3_ex7
grade_lab3_ex7(sampler_result, result_prob_dist) # Expected result type - sampler_result: SamplerResult, result_prob_dist: ProbDistribution
Submitting your answer. Please wait... Congratulations 🎉! Your answer is correct and has been submitted.
# Compute cost of the obtained result and display result
TSPCost(energy = energy, result_bitstring = result_bitstring[7:], adj_matrix = adj_matrix)
------------------- Energy: -4750.24525 Tsp objective: 130.75475000000006 Feasibility: True Solution Vector: [1, 0, 2] Solution Objective: 130.0 -------------------
# Plot convergence
PlotGraph(ideal = ideal, mean = [mean])
Next, let's try to reduce our search subspace even further! We shall now take into consideration the constraint specified in Equation (2). Here, unlike the first PQC, we need to add more correlations among the qubits since we will be mapping more variables across the whole matrix representation. Since we have variables that appear both in the first and second line, we can no longer realize the constraints by just tensor products. Hence we will entangle the parametrized $W$ gates using $CNOT$ gates to achieve the entanglement as mentioned. Note, the qubit mapping here is in the column-major ordered for the code and the columns corresponding to the problem formulation in our case, so we shall formulate according to this format. The illustrations in the paper may have row/column major interchanged so do keep an eye out for the notation followed!
We shall break down the protocol into a few steps and build up the routine for a 3 node graph. As mentioned above, we shall be taking into account not only the first line but also the second line of constraints to further reduce the subspace. To add in correlations among the qubits mapped to the variables (recall the mapping is column-major ordered), we will now create an entangled routine to complete the protocol.
Let us first create a PQC that satisfies both of the constraints $\sum_{p=1}^N x_{1,p}=1$ and $\sum_{v=1}^N x_{1,v}=1$. Recall, here $v$ refers to the vertex and $p$ represents its order in a prospective cycle. Building for a 3 node graph, let us first build the 'L' constraint circuit as shown in the paper!
from qiskit import QuantumRegister, ClassicalRegister
from qiskit.circuit import ParameterVector, QuantumCircuit
l_circuit = QuantumCircuit(n**2)
theta = ParameterVector('theta',length=(n-1)*n-1)
# L-Shaped Entangler
l_circuit.ry(theta[0],0)
l_circuit.cx(0,1)
l_circuit.cx(0,3)
l_circuit.barrier()
#W_phi_p
l_circuit.x(1)
l_circuit.ry(theta[1],2)
l_circuit.cz(1,2)
l_circuit.ry(-theta[1],2)
l_circuit.cx(2,1)
#W_phi_v
l_circuit.x(3)
l_circuit.ry(theta[2],6)
l_circuit.cz(3,6)
l_circuit.ry(-theta[2],6)
l_circuit.cx(6,3)
l_circuit.draw("mpl")
Next, now we need to encode the remaining constraints excluding the ones we applied. Since, the constraints for the remaining part are already determined, they can be read in a similar way as the previous problem by applying the corresponding $CNOT$ gates followed by the parametrized $W$ state gates on the qubits mapped variables as shown below!
from qiskit import QuantumRegister, ClassicalRegister
from qiskit.circuit import ParameterVector, QuantumCircuit
r_circuit = QuantumCircuit(n**2)
theta = ParameterVector('theta',length=(n-1)*n-1)
# Two blocks remaining for N = 3
# Block 1 ---------
r_circuit.x(4)
r_circuit.ry(theta[3],5)
r_circuit.cz(4,5)
r_circuit.ry(-theta[3],5)
r_circuit.cx(5,4)
# Block 2 ---------
r_circuit.x(7)
r_circuit.ry(theta[4],8)
r_circuit.cz(7,8)
r_circuit.ry(-theta[4],8)
r_circuit.cx(8,7)
r_circuit.draw("mpl")
One final step before we join the circuits is making the bridge between L-shape and original entangler. We shall construct it as follows:
from qiskit import QuantumRegister, ClassicalRegister
from qiskit.circuit import ParameterVector, QuantumCircuit
bridge_circuit = QuantumCircuit(n**2)
theta = ParameterVector('theta',length=(n-1)*n-1)
bridge_circuit.cx(3,4)
bridge_circuit.cx(6,7)
bridge_circuit.draw(output="mpl")
Now we proceed to join all the three circuits as shown in the figure below to complete the routine. We will run the constructed circuit using the RunVQE function which we created in the previous exercises to evaluate the result!
from qiskit import QuantumRegister, ClassicalRegister
from qiskit.circuit import ParameterVector, QuantumCircuit
model_2 = QuantumCircuit(n**2)
theta = ParameterVector('theta',length=(n-1)*n-1)
model_2 = l_circuit.compose(bridge_circuit).compose(r_circuit)
model_2.draw("mpl")
Calling the RunVQE function to run the circuit for our graph with the following model:
# %%time
# # # #This is not a graded exercise. Run this cell when you have spsare time.
# # Define optimizer
# optimizer = SPSA()
# # Define Initial point
# np.random.seed(10)
# init_2 = np.random.rand((n-1)*n-1) * 2 * np.pi
# # Call RunVQE
# with Session(service = service, backend = backend):
# result_m2, mean_m2 = RunVQE(Estimator(options=options), model_2, optimizer, qubitOp, init=init_2)
# print(result_m2, "\n")
# # #This is not a graded exercise. Run this cell when you have spsare time.
# optimal_m2 = result_m2.optimal_circuit.measure_all(inplace=False)
# %%time
# # # #This is not a graded exercise. Run this cell when you have spsare time.
# # Compute values
# with Session(service=service, backend=backend):
# energy_m2, sampler_result_m2, result_val_m2, result_m2_bitstring = BestBitstring(result=result_m2, optimal_circuit=optimal_m2)
# print("Optimal bitstring = ", result_m2_bitstring[7:])
# # # #This is not a graded exercise. Run this cell when you have spsare time.
# # Compute cost of the obtained result and display result
# TSPCost(energy=energy_m2, result_bitstring=result_m2_bitstring[7:], adj_matrix=adj_matrix)
# Plot convergence
#PlotGraph(mean = [mean,mean_m2], ideal = ideal)
We can see the faster convergence with the newer model! We shall now look at the fourth case extending on the implementation of the second case here. We shall be skipping the third case in the interest of keeping the lab length in check, but please feel free to check out the "parameter sharing" approach of constructing the PQC shared in the paper that aims to keep the circuit depth in check.
Now let us consider the final case for the problem. For the fourth PQC, we consider all constraints to completely exclude the infeasible answers in the above image. Thus, the set of the bases of the quantum states includes only feasible answers. We shall now work on an intuition as to how to exclude all the infeasible results. If we see the 2D grid above, based on the constraints we define, the feasible answers can be interpreted in the form of permutation matrices.
A permutation matrix is a sqaure matrix that is made up of only 1s and 0s and each row and column have exactly one 1 only. You may know the matrix $I$, which is called the identity matrix. It is a special case of permutation matrices.
Now that we know the definition of a permutation matrix, we can count the number of feasible answers that exist. If n = 2, there are 2 feasible solutions. Below you can check out what these matrices are.
Reflecting these two feasible answers, the circuit below is an example which satisfies all the constraints when n equals 2.
n = 2
perm_circuit = QuantumCircuit(n**2)
theta = ParameterVector('theta',length=(n-1)*n-1)
perm_circuit.x(0)
perm_circuit.ry(theta[0],1)
perm_circuit.cz(0,1)
perm_circuit.ry(-theta[0],1)
perm_circuit.cx(1,0)
perm_circuit.cx(1,2)
perm_circuit.cx(0,3)
perm_circuit.draw("mpl")
Such a PQC for arbitrary $N$ can be constructed in a recursive manner. Because we know the circuit when n = 2, we can construct circuits for the cases $n = 3, 4, ... k$. Let's make the $k \times k$ permutation matrices from the $(k-1) \times (k-1)$ ones.
As the steps above, we can create the desired quantum states as follows:
Let's start on the previous circuit while adding 5 more qubits. Note that CNOT gates' position is a little different because now $q_2$ in the circuit means $q_{3,1}$.
Exercise 8:
Build up the model_3 circuit for a 3 node graph problem to accommodate all the constraints of the problem as shown in the paper. We shall use the previous helper functions to run our problem.
n = 3
model_3 = QuantumCircuit(n**2)
theta = ParameterVector('theta',length=(n-1)*n//2)
model_3.x([0,6])
model_3.ry(theta[0],1)
model_3.ry(theta[1],7)
model_3.ry(theta[2],8)
model_3.cz(0,1)
model_3.cz(6,7)
model_3.ry(-theta[0],1)
model_3.ry(-theta[1],7)
model_3.cz(7,8)
model_3.ry(-theta[2],8)
model_3.cnot(1,0)
model_3.cnot(1,3)
model_3.cnot(0,4)
model_3.cnot(7,6)
model_3.cnot(8,7)
model_3.cswap(6,0,2)
model_3.cswap(6,3,5)
model_3.cswap(7,1,2)
model_3.cswap(7,4,5)
model_3.draw("mpl")
The circuit above represents exactly the superposition of bases of 6 feasible answers. Now we shall call our helper function RunVQE to execute this routine
%%time
# Define optimizer
optimizer = SPSA(maxiter=50)
# Define Initial point
np.random.seed(10)
init_3 = np.random.rand((n-1)*n//2) * 2 * np.pi
#### Enter your code below ####
# Call RunVQE
with Session(service = service, backend = backend):
result_m3, mean_m3 = RunVQE(estimator=Estimator(options=options),model=model_3,operator=qubitOp,optimizer=optimizer,init=init_3)
#### Enter your code above ####
CPU times: user 2.28 s, sys: 187 ms, total: 2.46 s Wall time: 5min 46s
optimal_m3 = result_m3.optimal_circuit.measure_all(inplace=False)
%%time
# Compute values
with Session(service=service, backend=backend):
energy_m3, sampler_result_m3, result_val_m3, result_m3_bitstring = BestBitstring(result=result_m3, optimal_circuit=optimal_m3)
print("Optimal bitstring = ", result_m3_bitstring[7:])
Optimal bitstring = [0 1 0 1 0 0 0 0 1] CPU times: user 46.9 ms, sys: 3.05 ms, total: 49.9 ms Wall time: 22.3 s
from qc_grader.challenges.fall_2022 import grade_lab3_ex8
grade_lab3_ex8(model_3, result_m3) # Expected result type - model_3: QuantumCircuit, result_m3: VQEResult
Submitting your answer. Please wait... Congratulations 🎉! Your answer is correct and has been submitted.
# Compute cost of the obtained result and display result
TSPCost(energy = energy_m3, result_bitstring = result_m3_bitstring[7:], adj_matrix = adj_matrix)
------------------- Energy: -4751.000000000001 Tsp objective: 129.9999999999991 Feasibility: True Solution Vector: [1, 0, 2] Solution Objective: 130.0 -------------------
# Plot convergence
PlotGraph(mean = [mean,mean_m3], ideal = ideal)
# Plot all three if you have run the L-shaped optional cell
# PlotGraph(mean = [mean,mean_m2,mean_m3], ideal = ideal)
You will see a flat lined result close to the ideal result with an ideal simulator because it contains essentially only the feasible solution subspace which converges very fast in an ideal scenario.
Remember our example for a 4 node graph which failed to converge? Now it's time to give it another go! Formulate a PQC for a 4 node graph in order to attempt a convergence for a given 4 node graph!
### DO NOT CHANGE BELOW
n = 4
num_qubits = n**2
# the target graph
tsp = Tsp.create_random_instance(n, seed=123)
colors = ["tab:red" for node in tsp.graph.nodes]
pos = [tsp.graph.nodes[node]["pos"] for node in tsp.graph.nodes]
draw_graph(tsp.graph, colors, pos)
### DO NOT CHANGE ABOVE
#### Enter your code below ####
qp = tsp.to_quadratic_program()
qp2qubo = QuadraticProgramToQubo()
qubo = qp2qubo.convert(qp)
qubitOp_4, _ = qubo.to_ising()
#### Enter your code above ####
model = tsp_entangler(4)
%%time
from qiskit.utils import algorithm_globals
algorithm_globals.random_seed = 1024
options = Options(simulator={"seed_simulator": 42},resilience_level=0)
# Define optimizer
optimizer = SPSA(maxiter=50)
# Define Initial point
np.random.seed(10)
init = np.random.rand((n-1)*n) * 2 * np.pi
# Define backend
backend = service.backends(simulator=True)[0]
# Call RunVQE
with Session(service = service, backend = backend):
result_4, mean = RunVQE(Estimator(options=options), model, optimizer, qubitOp_4, init=init)
# Print result
print(result_4)
--------------------------------------------------------------------------- MaxRetryError Traceback (most recent call last) File /opt/conda/lib/python3.8/site-packages/requests/adapters.py:440, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies) 439 if not chunked: --> 440 resp = conn.urlopen( 441 method=request.method, 442 url=url, 443 body=request.body, 444 headers=request.headers, 445 redirect=False, 446 assert_same_host=False, 447 preload_content=False, 448 decode_content=False, 449 retries=self.max_retries, 450 timeout=timeout 451 ) 453 # Send the request. 454 else: File /opt/conda/lib/python3.8/site-packages/urllib3/connectionpool.py:876, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, **response_kw) 875 log.debug("Retry: %s", url) --> 876 return self.urlopen( 877 method, 878 url, 879 body, 880 headers, 881 retries=retries, 882 redirect=redirect, 883 assert_same_host=assert_same_host, 884 timeout=timeout, 885 pool_timeout=pool_timeout, 886 release_conn=release_conn, 887 chunked=chunked, 888 body_pos=body_pos, 889 **response_kw 890 ) 892 return response File /opt/conda/lib/python3.8/site-packages/urllib3/connectionpool.py:876, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, **response_kw) 875 log.debug("Retry: %s", url) --> 876 return self.urlopen( 877 method, 878 url, 879 body, 880 headers, 881 retries=retries, 882 redirect=redirect, 883 assert_same_host=assert_same_host, 884 timeout=timeout, 885 pool_timeout=pool_timeout, 886 release_conn=release_conn, 887 chunked=chunked, 888 body_pos=body_pos, 889 **response_kw 890 ) 892 return response [... skipping similar frames: HTTPConnectionPool.urlopen at line 876 (2 times)] File /opt/conda/lib/python3.8/site-packages/urllib3/connectionpool.py:876, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, **response_kw) 875 log.debug("Retry: %s", url) --> 876 return self.urlopen( 877 method, 878 url, 879 body, 880 headers, 881 retries=retries, 882 redirect=redirect, 883 assert_same_host=assert_same_host, 884 timeout=timeout, 885 pool_timeout=pool_timeout, 886 release_conn=release_conn, 887 chunked=chunked, 888 body_pos=body_pos, 889 **response_kw 890 ) 892 return response File /opt/conda/lib/python3.8/site-packages/urllib3/connectionpool.py:866, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, **response_kw) 865 try: --> 866 retries = retries.increment(method, url, response=response, _pool=self) 867 except MaxRetryError: File /opt/conda/lib/python3.8/site-packages/qiskit_ibm_runtime/api/session.py:107, in PostForcelistRetry.increment(self, method, url, response, error, _pool, _stacktrace) 98 logger.debug( 99 "Retrying method=%s, url=%s, status=%s, error=%s, data=%s, headers=%s", 100 method, (...) 105 headers, 106 ) --> 107 return super().increment( 108 method=method, 109 url=url, 110 response=response, 111 error=error, 112 _pool=_pool, 113 _stacktrace=_stacktrace, 114 ) File /opt/conda/lib/python3.8/site-packages/urllib3/util/retry.py:592, in Retry.increment(self, method, url, response, error, _pool, _stacktrace) 591 if new_retry.is_exhausted(): --> 592 raise MaxRetryError(_pool, url, error or ResponseError(cause)) 594 log.debug("Incremented Retry for (url='%s'): %r", url, new_retry) MaxRetryError: HTTPSConnectionPool(host='runtime-us-east.quantum-computing.ibm.com', port=443): Max retries exceeded with url: /jobs/cdp9944hl0189c638gdg (Caused by ResponseError('too many 524 error responses')) During handling of the above exception, another exception occurred: RetryError Traceback (most recent call last) File /opt/conda/lib/python3.8/site-packages/qiskit_ibm_runtime/api/session.py:269, in RetrySession.request(self, method, url, bare, **kwargs) 268 self._log_request_info(final_url, method, kwargs) --> 269 response = super().request(method, final_url, headers=headers, **kwargs) 270 response.raise_for_status() File /opt/conda/lib/python3.8/site-packages/requests/sessions.py:529, in Session.request(self, method, url, params, data, headers, cookies, files, auth, timeout, allow_redirects, proxies, hooks, stream, verify, cert, json) 528 send_kwargs.update(settings) --> 529 resp = self.send(prep, **send_kwargs) 531 return resp File /opt/conda/lib/python3.8/site-packages/requests/sessions.py:645, in Session.send(self, request, **kwargs) 644 # Send the request --> 645 r = adapter.send(request, **kwargs) 647 # Total elapsed time of the request (approximately) File /opt/conda/lib/python3.8/site-packages/requests/adapters.py:510, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies) 509 if isinstance(e.reason, ResponseError): --> 510 raise RetryError(e, request=request) 512 if isinstance(e.reason, _ProxyError): RetryError: HTTPSConnectionPool(host='runtime-us-east.quantum-computing.ibm.com', port=443): Max retries exceeded with url: /jobs/cdp9944hl0189c638gdg (Caused by ResponseError('too many 524 error responses')) The above exception was the direct cause of the following exception: RequestsApiError Traceback (most recent call last) File /opt/conda/lib/python3.8/site-packages/qiskit/algorithms/minimum_eigensolvers/vqe.py:254, in VQE._get_evaluate_energy.<locals>.evaluate_energy(parameters) 253 job = self.estimator.run(batch_size * [ansatz], batch_size * [operator], parameters) --> 254 estimator_result = job.result() 255 except Exception as exc: File /opt/conda/lib/python3.8/site-packages/qiskit_ibm_runtime/runtime_job.py:197, in RuntimeJob.result(self, timeout, decoder) 196 if self._results is None or (_decoder != self._final_result_decoder): --> 197 self.wait_for_final_state(timeout=timeout) 198 if self._status == JobStatus.ERROR: File /opt/conda/lib/python3.8/site-packages/qiskit_ibm_runtime/runtime_job.py:277, in RuntimeJob.wait_for_final_state(self, timeout) 276 time.sleep(0.1) --> 277 status = self.status() 278 except futures.TimeoutError: File /opt/conda/lib/python3.8/site-packages/qiskit_ibm_runtime/runtime_job.py:233, in RuntimeJob.status(self) 228 """Return the status of the job. 229 230 Returns: 231 Status of this job. 232 """ --> 233 self._set_status_and_error_message() 234 return self._status File /opt/conda/lib/python3.8/site-packages/qiskit_ibm_runtime/runtime_job.py:374, in RuntimeJob._set_status_and_error_message(self) 373 if self._status not in JOB_FINAL_STATES: --> 374 response = self._api_client.job_get(job_id=self.job_id()) 375 self._set_status(response) File /opt/conda/lib/python3.8/site-packages/qiskit_ibm_runtime/api/clients/runtime.py:210, in RuntimeClient.job_get(self, job_id) 202 """Get job data. 203 204 Args: (...) 208 JSON response. 209 """ --> 210 response = self._api.program_job(job_id).get() 211 logger.debug("Runtime job get response: %s", response) File /opt/conda/lib/python3.8/site-packages/qiskit_ibm_runtime/api/rest/program_job.py:51, in ProgramJob.get(self) 46 """Return program job information. 47 48 Returns: 49 JSON response. 50 """ ---> 51 return self.session.get(self.get_url("self")).json() File /opt/conda/lib/python3.8/site-packages/requests/sessions.py:542, in Session.get(self, url, **kwargs) 541 kwargs.setdefault('allow_redirects', True) --> 542 return self.request('GET', url, **kwargs) File /opt/conda/lib/python3.8/site-packages/qiskit_ibm_runtime/api/session.py:292, in RetrySession.request(self, method, url, bare, **kwargs) 291 raise IBMNotAuthorizedError(message) from ex --> 292 raise RequestsApiError(message, status_code) from ex 294 return response RequestsApiError: "HTTPSConnectionPool(host='runtime-us-east.quantum-computing.ibm.com', port=443): Max retries exceeded with url: /jobs/cdp9944hl0189c638gdg (Caused by ResponseError('too many 524 error responses'))" The above exception was the direct cause of the following exception: AlgorithmError Traceback (most recent call last) File <timed exec>:18, in <module> Input In [75], in RunVQE(estimator, model, optimizer, operator, init) 15 vqe = VQE(estimator=estimator,ansatz=model,optimizer=optimizer,initial_point=init) 16 # 17 # 18 # build your code here (...) 21 22 # Compute minimum_eigenvalue ---> 23 result = vqe.compute_minimum_eigenvalue(operator) 25 #### Enter your code above #### 27 return result, history["mean"] File /opt/conda/lib/python3.8/site-packages/qiskit/algorithms/minimum_eigensolvers/vqe.py:191, in VQE.compute_minimum_eigenvalue(self, operator, aux_operators) 186 else: 187 # we always want to submit as many estimations per job as possible for minimal 188 # overhead on the hardware 189 was_updated = _set_default_batchsize(self.optimizer) --> 191 optimizer_result = self.optimizer.minimize( 192 fun=evaluate_energy, x0=initial_point, jac=evaluate_gradient, bounds=bounds 193 ) 195 # reset to original value 196 if was_updated: File /opt/conda/lib/python3.8/site-packages/qiskit/algorithms/optimizers/spsa.py:556, in SPSA.minimize(self, fun, x0, jac, bounds) 554 iteration_start = time() 555 # compute update --> 556 fx_estimate, update = self._compute_update(fun, x, k, next(eps), lse_solver) 558 # trust region 559 if self.trust_region: File /opt/conda/lib/python3.8/site-packages/qiskit/algorithms/optimizers/spsa.py:487, in SPSA._compute_update(self, loss, x, k, eps, lse_solver) 484 num_samples = self.resamplings 486 # accumulate the number of samples --> 487 value, gradient, hessian = self._point_estimate(loss, x, eps, num_samples) 489 # precondition gradient with inverse Hessian, if specified 490 if self.second_order: File /opt/conda/lib/python3.8/site-packages/qiskit/algorithms/optimizers/spsa.py:464, in SPSA._point_estimate(self, loss, x, eps, num_samples) 461 delta1 = deltas1[i] 462 delta2 = deltas2[i] if self.second_order else None --> 464 value_sample, gradient_sample, hessian_sample = self._point_sample( 465 loss, x, eps, delta1, delta2 466 ) 467 value_estimate += value_sample 468 gradient_estimate += gradient_sample File /opt/conda/lib/python3.8/site-packages/qiskit/algorithms/optimizers/spsa.py:425, in SPSA._point_sample(self, loss, x, eps, delta1, delta2) 422 self._nfev += 2 424 # batch evaluate the points (if possible) --> 425 values = _batch_evaluate(loss, points, self._max_evals_grouped) 427 plus = values[0] 428 minus = values[1] File /opt/conda/lib/python3.8/site-packages/qiskit/algorithms/optimizers/spsa.py:744, in _batch_evaluate(function, points, max_evals_grouped, unpack_points) 742 results += _as_list(function(*batch)) 743 else: --> 744 results += _as_list(function(batch)) 746 return results File /opt/conda/lib/python3.8/site-packages/qiskit/algorithms/minimum_eigensolvers/vqe.py:256, in VQE._get_evaluate_energy.<locals>.evaluate_energy(parameters) 254 estimator_result = job.result() 255 except Exception as exc: --> 256 raise AlgorithmError("The primitive job to evaluate the energy failed!") from exc 258 values = estimator_result.values 260 if self.callback is not None: AlgorithmError: 'The primitive job to evaluate the energy failed!'
%%time
from qiskit.utils import algorithm_globals
algorithm_globals.random_seed = 1024
options = Options(simulator={"seed_simulator": 42},resilience_level=0)
# Define optimizer
optimizer = SPSA(maxiter=100)
# Define Initial point
np.random.seed(10)
init = np.random.rand((n-1)*n) * 2 * np.pi
# Define backend
backend = service.backends(simulator=True)[0]
# Call RunVQE
with Session(service = service, backend = backend):
result_4, mean = RunVQE(Estimator(), model, optimizer, qubitOp_4, init=init)
# Print result
print(result_4)
{ 'aux_operators_evaluated': None,
'cost_function_evals': 200,
'eigenvalue': -51450.829,
'optimal_circuit': <qiskit.circuit.quantumcircuit.QuantumCircuit object at 0x7fac3f01ebe0>,
'optimal_parameters': { ParameterVectorElement(theta[8]): 4.000723992577986,
ParameterVectorElement(theta[9]): -0.0015268286153453173,
ParameterVectorElement(theta[6]): 1.5800022432326093,
ParameterVectorElement(theta[3]): 1.574533018138417,
ParameterVectorElement(theta[0]): 4.697714469145179,
ParameterVectorElement(theta[5]): 1.5828090187805681,
ParameterVectorElement(theta[1]): 1.5702229926394178,
ParameterVectorElement(theta[2]): 3.142003904248961,
ParameterVectorElement(theta[7]): 6.270739390871984,
ParameterVectorElement(theta[11]): 3.736549300151775,
ParameterVectorElement(theta[4]): 4.702630301223836,
ParameterVectorElement(theta[10]): 4.107171184051593},
'optimal_point': array([ 4.69771447e+00, 1.57022299e+00, 3.14200390e+00, 1.57453302e+00,
4.70263030e+00, 1.58280902e+00, 1.58000224e+00, 6.27073939e+00,
4.00072399e+00, -1.52682862e-03, 4.10717118e+00, 3.73654930e+00]),
'optimal_value': -51450.829,
'optimizer_evals': None,
'optimizer_result': <qiskit.algorithms.optimizers.optimizer.OptimizerResult object at 0x7fabf04b0ee0>,
'optimizer_time': 930.4119226932526}
CPU times: user 5.81 s, sys: 390 ms, total: 6.2 s
Wall time: 15min 30s
optimal = result_4.optimal_circuit.measure_all(inplace=False)
%%time
# Compute values
with Session(service=service, backend=backend):
energy, sampler_result, result_val, bitstring = BestBitstring(result=result_4, optimal_circuit=optimal)
print("Optimal bitstring = ", bitstring)
Optimal bitstring = [0 0 0 1 0 0 1 0 1 0 0 0 0 1 0 0] CPU times: user 60.7 ms, sys: 3.59 ms, total: 64.3 ms Wall time: 1.18 s
# Plot convergence
ee = NumPyMinimumEigensolver()
result = ee.compute_minimum_eigenvalue(qubitOp_4)
energy_numpy = result.eigenvalue.real
PlotGraph(mean = [mean], ideal = energy_numpy)
from qc_grader.challenges.fall_2022 import grade_lab3_ex9
grade_lab3_ex9(result_4, bitstring)
Submitting your answer. Please wait...
Congratulations 🎉! Your answer is correct.
You've successfully carved a clear path for your slingshot by targeting and collecting
certain orbiting objects with your drone.
You are almost ready to perform your slingshot. As you move your starship into position,
you wonder curiously about the other you out there.
Have they ceased to exist now that you haven't followed their path by selecting this planet?
Do they somehow exist in an alternate reality?
You can never know.
Nonetheless, you look out the window at the planets you didn't select, and whisper an aching,
heartfelt, "Thank you."
We shall now explore the effect of Zero Noise Extrapolation (ZNE) on our routine. Recall in Lab 1 we explored that error mitigation strategies on runtime can be switched on when specifying the resilience level and higher levels generate more accurate results, at the expense of longer processing times.. Here we shall now specify a noise model and resilience levels to play around with ZNE.
Please note: This feature is still in beta and results might not be as expected. The module is in continuous development scheduled for stable release in the future so feel free to experiment around with this routine!
We shall use the FakeGuadalupe backend to simulate our noisy system. Setting the options for a baseline run:
from qiskit.providers.fake_provider import FakeGuadalupe
from qiskit_aer.noise import NoiseModel
fake_backend = FakeGuadalupe()
noise_model = NoiseModel.from_backend(fake_backend)
options = Options(
simulator={
"noise_model": noise_model,
"seed_simulator": 42,
},
resilience_level=0
)
# Define 3 Node graph for our learning exercises
n = 3
num_qubits = n**2
tsp = Tsp.create_random_instance(n, seed=250)
adj_matrix = nx.to_numpy_matrix(tsp.graph)
colors = ["tab:red" for node in tsp.graph.nodes]
pos = [tsp.graph.nodes[node]["pos"] for node in tsp.graph.nodes]
print("Learning problem 3 Node graph: ")
draw_graph(tsp.graph, colors, pos)
# Define quadratic program for 4 node graph
qp = tsp.to_quadratic_program()
# Define Ising operator for 3 Node graph
qp2qubo = QuadraticProgramToQubo()
qubo = qp2qubo.convert(qp)
qubitOp, offset = qubo.to_ising()
# Making the Hamiltonian in its full form and getting the lowest eigenvalue and eigenvector
ee = NumPyMinimumEigensolver()
result = ee.compute_minimum_eigenvalue(qubitOp)
ideal = result.eigenvalue.real
print(f"Ideal energy: {energy_numpy}")
Learning problem 3 Node graph: Ideal energy: -51520.0
Now to get a baseline result, we shall run all our models we made above with a noisy backend. Here in the run options above, we shall specify the noise model and noise parameters to be passed on to our simulator backend. We shall be using 'FakeGuadalupe' as our backend.
Do note, this can take well over 15 minutes to run so feel free to take some rest and stretch your legs as this is being processed!
%%time
optimizer = SPSA(maxiter= 50)
init = np.random.rand((n-1)*n) * 2 * np.pi
model_1 = tsp_entangler(3)
with Session(service = service, backend = backend):
print("Running model 1......")
result_n, mean_n = RunVQE(Estimator(options=options), model_1, optimizer, qubitOp, init=init)
print("Model 1 simulation complete")
Running model 1...... Model 1 simulation complete CPU times: user 29.8 s, sys: 343 ms, total: 30.1 s Wall time: 11min 29s
%%time
with Session(service = service, backend = backend):
print("Running model 2......")
result_m2_n, mean_m2_n = RunVQE(Estimator(options=options), model_2, optimizer, qubitOp, init=np.random.rand((n-1)*n-1) * 2 * np.pi)
print("Model 2 simulation complete")
Running model 2...... Model 2 simulation complete CPU times: user 30.1 s, sys: 310 ms, total: 30.4 s Wall time: 30min 27s
%%time
n=3
init_3 = np.random.rand((n-1)*n//2) * 2 * np.pi
with Session(service = service, backend = backend):
print("Running model 3......")
result_m3_n, mean_m3_n = RunVQE(Estimator(options=options), model_3, optimizer, qubitOp, init=init_3)
print("Model 3 simulation complete")
Running model 3...... Model 3 simulation complete CPU times: user 30.1 s, sys: 332 ms, total: 30.4 s Wall time: 15min 14s
# Plot convergence
PlotGraph(mean = [mean_n,mean_m2_n,mean_m3_n], ideal = ideal)
Notice the result graphs are no longer converging to the ideal point as before due to the noise model applied. Now let's try to get it as close to the ideal value as possible by using Digital ZNE on Qiskit Runtime!
(Adapted from: Configure error mitigation release notes)
The Estimator interface lets users seamlessly work with the variety of error mitigation methods to reduce error in expectation values of observables. Below is an example of leveraging Zero Noise Extrapolation by simply setting resilience_level=2. The default settings sample the expectation value at three noise factors, leading to a roughly 3x overhead when employing this resilience level.
You can tune advanced options to configure your resilience strategy further. These methods can be used alongside resilience levels where you change the specific options of interest and let your previously set resilience level manage the rest. Please note, this is a beta release and as a part of the beta release of the resilience options, users will be able configure ZNE by using the following advanced options below:
| Options | Inputs | Desctiption |
|---|---|---|
| options.resilience.noise_amplifier(Optional[str]): Amplification Strategy | TwoQubitAmplifier [Default] |
Amplifies noise of all two qubit gates by performing local gate folding. |
CxAmplifier |
Amplifies noise of all CNOT gates by performing local gate folding. | |
LocalFoldingAmplifer |
Amplifies noise of all gates by performing local gate folding. | |
GlobalFoldingAmplifier |
Amplifies noise of the input circuit by performing global folding of the entire input circuit. | |
| options.resilience.noise_factors((Optional[Sequence[float]]) | (1, 3, 5) [Default] |
Noise amplification factors, where 1 represents the baseline noise. They all need to be greater than or equal to the baseline. |
| options.resilience.extrapolator(Optional[str]) | LinearExtrapolator [Default] |
Polynomial extrapolation of degree one. |
QuadraticExtrapolator |
Polynomial extrapolation of degree two and lower. | |
CubicExtrapolator |
Polynomial extrapolation of degree three and lower. | |
QuarticExtrapolator |
Polynomial extrapolation of degree four and lower. |
# Set options for Digital ZNE routine
options = Options(
simulator={
"noise_model": noise_model,
"seed_simulator": 42,
},
)
# Settings to set noise_factors, noise_amplifier and extrapolator
options.resilience.noise_factors = (1, 2, 3, 4) # Set noise_factors
options.resilience.noise_amplifer = 'CxAmplifer' # Set noise_amplifiers
options.resilience.extrapolator = 'QuadraticExtrapolator' # Set noise_extrapolators
%%time
optimizer = SPSA(maxiter=10)
with Session(service = service, backend = backend):
print("Running model 1......")
result_zn, mean_zn = RunVQE(Estimator(options=options), model_1, optimizer, qubitOp, init=init)
print("Model 1 simulation complete")
print("Running model 2......")
result_m2_zn, mean_m2_zn = RunVQE(Estimator(options=options), model_2, optimizer, qubitOp, init=np.random.rand((n-1)*n-1) * 2 * np.pi)
print("Model 2 simulation complete")
print("Running model 3......")
result_m3_zn, mean_m3_zn = RunVQE(Estimator(options=options), model_3, optimizer, qubitOp, init=init_3)
print("Model 3 simulation complete")
Running model 1...... Model 1 simulation complete Running model 2...... Model 2 simulation complete Running model 3...... Model 3 simulation complete CPU times: user 28.1 s, sys: 119 ms, total: 28.2 s Wall time: 21min 26s
# Plot convergence
PlotGraph(mean = [mean_zn,mean_m2_zn,mean_m3_zn], ideal = energy_numpy)
Created by: Astri Cornish, Dayeong Kang, Desiree Vogt-Lee, Vishal Bajpe, Yuri Kobayashi
Advisor: Takashi Imamichi, Atsushi Matsuo, Yudai Suzuki
Creative assets by: Radha Pyari Sandhir
Version: 1.0
from qiskit.tools.jupyter import *
%qiskit_version_table